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We study the Bose-Einstein condensation of an interacting gas with attractive interaction confined in 
a harmonic trap using a semiclassical two-fluid mean-field model. The condensed state is described 
by converged numerical solution of the Gross-Pitaevskii equation. By solving the system of coupled 
equations of this model iteratively we obtain converged results for the temperature dependencies of 
the condensate fraction, chemical potential, and internal energy for the Bose-Einstein condensate of 
7 Li atoms. 
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I. INTRODUCTION 

There has been recent experimental observation of 
Bose-Einstein (BE) condensation in dilute interacting 
bosonic atoms of 87 Rb @|§, 23 Na §, 7 Li [g,|, and : H 
H employing magnetic traps at ultra-low temperatures. 
The interaction among the atoms could be either attrac- 
tive or repulsive. Although, in both cases the condensate 
is well-described by the Gross-Pitaevskii (GP) equation 
(@,|), the nature of BE condensates in these two cases 
is of entirely different nature. In the repulsive case the 
number of atoms in the condensate can grow without 
bound, whereas this number is limited by a upper bound 
in the case of attractive interaction [p|-pd|. Of the ex- 
perimentally observed cases of BE condensation one has 
repulsion for 1 H, 23 Na, and 87 Rb atoms and attraction 
for 7 Li atoms. The existence of a maximum number of 
condensed atoms in the case of 7 Li has been noted ex- 
perimentally and is consistent with the prediction of the- 
oretical analysis based on the GP equation [|]j9] |l3| . In 
the attractive case the GP equation has no solution for 
the number of atoms larger than a critical number. 

The GP equation is a nonlinear Schrodinger equation 
and it is tedious to find its converged numerical solution 
Jl~2| [l7] . For the repulsive case, an approximate solution 
scheme of this equation, such as the one based on the 
Thomas-Fermi approximation, has frequently been used 
for a qualitative description of the condensate |l8| . How- 
ever, no such approximation scheme is known for the at- 
tractive case which requires an exact numerical solution 
of the GP equation. This makes the theoretical study of 
this case a more challenging task. 

There have been several comprehensive studies on the 
temperature dependencies of the thermodynamic observ- 
ables in the case of repulsive interaction using mean- 
field two-fluid models Jl8ppl|. One such mean- field 
scheme is provided by the so-called Popov approxima- 
tion |23] and has been considered by several authors 
p9|- pl| , p3 24 1 . The physical ingredients of these mean- 
field models JIJ|-|2l],^3],|4| are quite similar and they lead 
to similar numerical results in the case of weakly repul- 



sive interatomic interactions. 

We shall use a two-fluid mean-field model |l8| to study 
the temperature dependencies of the thermodynamic ob- 
servables of the condensate. For a condensate composed 
of 40000 trapped 87 Rb atoms the perturbative solution of 
the system of equations of this model converged rapidly 
and provided a satisfactory account of the condensate 
fraction, internal energy, and specific heat in agreement 
with experiment [QES]. It was also found that the lowest 
order solution already provided a very good approxima- 
tion. Later the same model has been used in one and two 
space dimensions 25 2q|. 

The above mean-field two-fluid model |1| is used in 
this work for a theoretical description of the BE conden- 
sation thermodynamics in the case of attractive interac- 
tion appropriate for 7 Li. Depending on the strength of 
the attractive potential, the condensate in the attractive 
case may consist of a few thousand atoms confined by the 
trap potential. For a fixed trap, the maximum number 
of atoms in the BE condensate with attractive interac- 
tion is inversely proportional to \a\ ]23| ], where a is the 
scattering length of two atoms. As the temperature is 
lowered below the critical temperature To of BE conden- 
sation, the condensate starts to form and finally at K 
all the available atoms (limited by the maximum num- 
ber mentioned above) form the condensate in the present 
model. 

The condensate wave function in the present model is 
described by the GP equation. In the repulsive case usu- 
ally some approximate solutions of the GP equation are 
used | jl8| , ^9p^ ]. Although we shall be using the itera- 
tive solution of the system of equations of the mean-field 
model, we shall employ a converged numerical solution 
of the GP equation in the present attractive case. As 
the GP equation is a nonlinear one, this amounts to a 
nontrivial modification of the calculational scheme. 

The plan of the paper is as follows. In Sec. II we 
present the mean-field two-fluid model. In Sec. Ill wc 
discuss the numerical scheme for its solution and present 
numerical results for 7 Li. Finally, in Sec. IV we present 
some concluding remarks. 
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II. MEAN-FIELD TWO-FLUID MODEL 

We consider a system of N bosons with attractive in- 
teraction at temperature T under the influence of a trap 
potential. The condensate is described by the following 
GP equa tion for the wave function ^(r) with eigenvalue 



^v 2 

2m 



f(r) = 0. 
(2.1) 



Here V ex ^(r) = mui 2 r 2 /2 is the spherically symmetric 
harmonic-oscillator trap potential, g = -iTrh 2 a/m the 
strength of the interatomic interaction, m the mass of 
a single atom, lu the angular frequency, a the atom-atom 
scattering length, p = p — po is t he c hemical potential, 
where po is the eigenvalue of Eq. (2J) for the harmonic 
oscillator potential alone in the absence of interatomic 
interaction (g = 0), and ni(r) represents the distribu- 
tion function of the noncondensed bosons. Although in 
actual experiment j^] the harmonic oscillator trap is not 
quite symmetric, the deviation from spherical symmetry 
is quite small. However, the converged numerical solu- 
tion of the GP equation (2.1) for a nonsymmetric trap is 
quite complicated numerically and hence for a qualitative 
description we consider a spherically symmetric trap in 
the present study. An attractive (repulsive) interaction 
corresponds to negative (positive) values of a and g. 

The noncondensed particles are treated as non- 
interacting bosons in an effective potential V e ^{r) = 

v ext( r ) + 2 .9 n i ( r ) + 2 3* 2 ( r ) @- Thermal averages are 
calculated with a standard Bose distribution of the non- 
condensed particles in chemical equilibrium with the con- 
densate governed by the same chemical potential p. In 
particular the density n\(r) is given by |27|] 



ni(r) 



1 



d 3 p 



(2ttH) 3 J exp[{p 2 /2m + V eff (r 



fi}/k B T]-r 
(2.2) 



wh ere k B is the Boltzmann constant. Equations (2.1) 
— (2.2) above are the principal equations of the present 
model. 

The total number of particles N of the system is given 

by 



N = N + 



p(E)dE 



exp[(£- p)/k B T] - 1' 



(2.3) 



where Nq = J ^> 2 (r)d 3 r is the total number of particles in 
the con den sate. The semiclassical density of states p(E) 
of Eq. (|U|) is given by || 



P(E) 



(2m) 



3/2 



An2?l3 Jv e S {r)<E 



E - V cS (r)d 3 



(2.4) 



The cri tica l temperature To is obtained as the solution 
of Eq. (2J3) with N Q and p set equal to 0. 

The average single-particle energy of the noncondensed 
particles is given by Q 



(£)nc = J 



Ep(E)dE 



exp[(E-n)/k B T\-l' 



(2.5) 



The kinetic energy of the condensate is assumed to be 
negligible and its interaction energy per particle is given 
by (E)c = (g/2) J 1 $> i (r)d 3 r. The quantity of experimen- 
tal interest is the average total energy (E) — [(E) n c(N — 
N )/2 + (E) c ]/N, which we calculate. 

The coupling of the nonlinear equation (^[l]) with other 
equations make the solution algorithm complicated and 
it is convenient to express this system of equation in di- 
mensionless units. This has advantage in the numeri- 
cal solution 0,0. We express energy in units of hoj, 
and length in units of the harmonic oscillator length 
a no = yjTi/(moj). We shall consider in this paper only 
the spherically symmetric ground state solution of the 
cond ensate with \&(r) = \&(r). Then the GP equation 
(|2.lD become 



dx 2 



Airrjni (x) ± 



2cj> 2 {x) 



2a 



(x)=0, (2.6) 



where rj = 4a/a n0 is the new dimensionless strength, 

x = r/a ho , *(r) = <j>(x) / \x J V^a^) , and a = fi/(huj). 

The positive (negative) sign in Eq. ( |2.6| ) corresponds 
to repulsive (attractive) interaction. The dimensionless 
density ni(x) = a^ Q ni(r) is defined by 



ni(x) 



k 2 dk 



2ir 2 J exp[{fc 2 /2 + ^ e ff(- T ) - - 1 



. (2-7) 



where k = pa^/fr, V c ff{ x ) — x 2 /2 + 2irr]ni(x) + 
2(/) 2 {x)/x 2 , t = k B T/(huj), a = p/(huj). For the spher- 
ically symmetric ground state a = a — q.q, where «o = 
po/(huj) = 1.5 is the energy eigenvalue of the harmonic 
oscillator potential alone (zero-point energy in units of 
hcu) in the absence of interatomic interaction. In the 
present consideration of the chemical potential finite-size 
effects are excluded jl^jljj. U sing these dimensionless 
variables the number equation (2.3) becomes 



N = N + 



p(e)de 



exp[(e — a)/t] — 1 



(2.8) 



where N Q = (4/|r?|) J (j) 2 (x)dx, e = E/(hu)), and the di- 
mensionless density p(e) = huip(E) is given by 



P(e) 



2y/2 



V eS ( X )<e 



^e-V eS {x)x 2 dx. (2.9) 



The above set of equations fl2.6| ) — ( |2.9[ ) are solved iter- 
atively using the converged numerical solution of the GP 
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equation (2.6). The iteration is started with n-y(x) = at 
a definite temperature w ith a trial value for the chemical 
potential a. Then Eq. ( |2.6D is solved and with its solu- 
tion the functions V g q(x) and fii( x) a re calculated. Using 
these new V g q(x) and n\(x) Eq. (2.6) is solved again and 
n\{x) and ^(a;) are recalculated. This iterative scheme 
is continued until final convergence is achieved. In each 
order of iteration we calculate the condensate fraction 
Nq/N and energy (E), in addition to the chemical po- 
tential a. This scheme is repeated until convergence is 
achieved. It is then verified if the number equation ( |2.S| ) 
is satisfied with this solution. If not, a new trial value 
for the chemical potential is employed. Once the number 
equation is satisfied the desired solution is obtained. 

Ne xt w e present the solution procedure of the GP equa- 
tion (2.6). The solution of this equation satisfies the fol- 



lowing boundary conditions 
0(0) = 0, 
lim </>(x) = Nc exp 
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4>'(0) = constant 
' x 2 , 1 



ln.x 



(2.10) 
(2.11) 



where Nc is a normalization constant. The deri vative 
of the wave function can be obtained from Eq. (2.11) 



and one obtains the following log-derivative of the wave 
function in the asymptotic region 



lim 



4>{x) 









-x + 




±1 






X 



(2.12) 



Equation ( |2.6| ) is integrated numerically for a given a 
by the four-point Runge-Kutta rule starting at th e ori - 
gin (x — 0) with the initial boundary condition (2.1C) 
with a trial </>'(0) in steps of dx = 0.001. Using Eq. 
(2.6) the integration is pr opaga ted to x = imax, where 
the asymptotic condition ( 2.12 ) is valid. The agreement 
between the numerically calculated log- deriv ative of the 
wave function and the theoretical result ( 2.12 ) is enforced 
to four significant figures. The maximum value of x up 
to which we need to integrate to obtain this precision is 
Xmax = 3. If for a trial <^'(0), this precision can not be 
obtained, a new value of (j)'(0) is to be chosen. The proce- 
dure is repeated until the converged solution is obtained. 



III. NUMERICAL RESULTS 

In our numerical study we would be interested only 
in the case of attractive interaction (negative rj). We 
consider the experimentally relevant situation of 7 Li . 
In this case the trap frequencies in the experiment of Rcf. 
pj along the X, Y, and Z directions are 150.6 Hz, 152.6 
Hz, and 131.5 Hz which lead to a maximum of about 1400 
trapped atoms. The deviation from spherical symmetry 
in this case is small and in order to have a qualitative 
understanding of the condensate we consider the trap 
to be spherically symmetric with -/Vniax = 1400. Wc 



reconfirm from the numerical solution of the GP equation 
that l2l 



max 



= 0.575. 



(3.1) 



Considering the known result for 7 Li, that iVmax = 1400 
H, we obtain from Eq. (B.l), 77 = —0.00164 and wc 
use this value of 77 in our numerical calculation. In 
our calculation we use three values of N, e.g., N = 
1300,1000, and 500 (N < iV max = 1400). First we 
calculate the critical temperatures in these cases from 
the number equation ( |2.8 ) and obtain the values To = 
l0.27?uj/kB,9Allicj/kB,7A7huj/kB for N = 1300,1000 
and 500, respectively. 

In the case of 7 Li, the estimated value of r}{= 
—0.00164) is quite small. This corresponds to a very weak 
coupling and we find that the lowest-order solution is 
graphically almost indistinguishable from the converged 
solution for the condensate fraction, chemical potential, 
and total energy. The estimated coupling 77 for 87 Rb 
is 0.0248 Ji],p|Jlq], which is more than ten times larger in 
magnitude than the coupling for 7 Li. In the case of 87 Rb, 
already the lowest order result was very good fll8| , f23"| . 
Hence the very rapid convergence of the present results 
is not quite unexpected. In the present study we only 
exhibit the converged result after two iterations. 
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Fig. 1. Converged condensate fraction Nq/N as a 
function of T/T for 77 = -0.00164 and N = 1300 (full 
line) and 500 (dashed-dotted line); for 77 = —0.01 and 
N = 200 (dashed-double-dotted line); and for trapped 
ideal Bose gas (dotted line). 

In Fig. 1 we present the temperature (T/Tq) de- 
pendence of the condensate fraction No/N for 7 Li with 
rj = -0.00164 for N = 1300 and 500. The result for 
N = 1000 is indistinguishable from that of N = 1300. 
The result for 77 = corresponding to an ideal Bose gas 
in a harmonic trap is also shown in this figure. For com- 
parison, the result in the case of a stronger attractive 
interaction for 77 = —0.01 and N — 200 is also shown. In 
the case of repulsive interaction, the result for N/Nq at 
a particular temperature T/Tq is smaller than that for 
an ideal Bose gas p3| . In the present case of attractive 
interaction, the result for N/Nq at a particular tempera- 
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ture T/Tq is larger than that for the ideal Bose gas. We 
also find from this figure that, as expected, the result for 
the stronger attractive interaction (77 = —0.01) deviates 
more from the ideal gas result than in the case of 7 Li 
(77 = —0.00164). We did not consider a mu ch stronger 
attractive coupling, as because of Eq. (3.1) this would 
correspond to an unacceptably small value for the num- 
ber of particles in the condensate. This number is already 
small for the case 77 — —0.01 considere d in this work. 

As the solution of the GP equation ( |2.6| ) in this case is 
nontrivial, we show in Fig. 2 the wave functions <j){x)/x 
for N = 1000 at temperatures T/T = 0, 0.4, 0.6, 0.8, 
and 0.9. As temperature decreases, the wave function 
is more pronounced corresponding to an increase in the 
number of particles Nq in the condensate given by the 
normalization Ao = (4/|r;|) J 4> 2 (x)dx. For other values 
of A, the wave functions are similar to those in Fig. 2 
and we do not show these wave functions here. 
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Figure 2 
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Fig. 2. GP wave function 4>{x)/x of Eq. (|2.6D for A - 
1000 at different temperatures^To for V = -0.00164. 
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Fig. 3. Chemical potential /i/(fcsTb) as a function 
of T/T for r] = -0.00164 and N = 1300 (full line), 
1000 (dashed line) and 500 (dashed-dotted line); and for 
r] = -0.01 and N = 200 (dashed-double-dotted line). 

In Fig. 3 we show the chemical potential of the sys- 
tem at different temperatures T/T Q for A = 1300, 1000, 
and 500 and 77 = —0.00164. For comparison we also 
show the result for the stronger attractive interaction 
with A = 200 and r\ = —0.01. Above the critical tem- 



perature T > Tq, the chemical potential for a trapped 
ideal Bose gas is negative and it becomes zero at the 
critical temperature Q . The same is true in the present 
simplified model where the very weakly interacting non- 
condensed gas is taken to be noninteracting. The only 
effect of interaction is considered via the condensate. The 
effect of the interaction among the atoms of the noncon- 
densed gas could be important for T just below To when 
there will be a large fraction of noncondensed gas and a 
small fraction of condensed gas. For T close to when 
most atoms are condensed such effect could be neglected. 
At present a correct description of BEC thermodynamics 
including the interaction among noncondensed atoms is 
beyond the scope of the present work. Hence, although 
it would be more appropriate to treat the noncondensed 
gas to be interacting, the effect of interaction in the treat- 
ment of the noncondensed gas is expected to be negligi- 
bly small in the present study of very weakly interacting 
Bose gas (77 = —0.0164) and is neglected. In the case 
of the ideal Bose gas, below the critical temperature the 
chemical potential is identically equal to zero. For the 
present case of the trapped Bose gas with attractive in- 
teraction, the chemical potential, after becoming zero at 
T = Tq from a negative value, becomes negative again 
as the temperature is reduced below the critical temper- 
ature. The chemical potential for the stronger attractive 
interaction (r/ = —0.01) deviates more from the trapped 
ideal gas result (jx = 0) below the critical temperature, 
than in the case of the weak attractive interaction of 7 Li 
(i] = —0.00164). In the case of repulsive interaction, the 
chemical potential becomes positive for temperatures be- 
low T = T @. 

Figure 4 




Fig. 4. Energy (E) /(Nk B T ) as a function of T/T for 
r) = -0.00164 and N = 1300 (full line) and 500 (dashed- 
dotted line); for 77 = 0.00164 and N = 1300 (dotted 
line); and the classical Maxwell-Boltzmann distribution 
(straight line). 

In Fig. 4 we plot the temperature dependence of en- 
ergy (E)/(Nk B T) for 77 = -0.00164 and N = 1300 and 
500. These two cases lead to almost identical energies 
as can be seen from Fig. 4. For comparison, in this 
figure we also show the result for the repulsive interac- 
tion for 77 = 0.00164 and A = 1300. The energy for 
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the trapped ideal Bose gas should lie between the en- 
ergies for the attractive and repulsive cases mentioned 
above. The classical Maxwell-Boltzmann result is also 
shown. For Bose-Einstein condensation to materialize 
the energy of the system should be lower that the clas- 
sical result. The energy in the attractive case is smaller 
than the cor resp onding repulsive case below critical tem- 
perature BPf. 



IV. CONCLUSION 

In conclusion, we studied the temperature dependence 
of condensate fraction, chemical potential, and total en- 
ergy for a trapped 7 Li gas consisting of 1300, 1000, and 
500 atoms with attractive interaction using a mean-field 
two-fluid model. The maximum number of atoms allowed 
in this case is 1400 H,|). We employed an iterative solu- 
tion scheme of the system of equations of this model as 
m Refs. flfQ. In the case of 7 Li the attractive interac- 
tion is very weak and the the system of equations leads 
to rapid convergence. The condensate was described by 
the converged numerical solution of the Gross-Pitaevskii 
equation |p^-p7|. As the interaction is weak, the results 
for condensate fraction, chemical potential, and energy 
of 7 Li are very close to the corresponding results for the 
trapped ideal Bose gas. However, the deviation from 
the result of the ideal Bose gas in case of 7 Li is in the 
opposite direction compared to the corresponding devia- 
tion in case of repulsive interaction [fl8[|23"|| . For example, 
in all cases the chemical potential is zero at the critical 
temperature. For ideal Bose gas it continues to be zero 
below the critical temperature. Below the critical tem- 
perature the chemical potential turns negative in case of 
7 Li, whereas it turns positive in case of 87 Rb |2j| where 
the interaction is repulsive. Although, in the present at- 
tractive case the number of particles in the condensate is 
small, the present results for the thermodynamic observ- 
ables are quite reasonable physically. This demonstates 
that the mean-field two-fluid models |[^,[L9 22 2^ used to 
study the thermodynamic observables for the BE conden- 
sate in the repulsive case are quite useful in the attractive 
case also. 
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